load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/home/Earth/rwhite//scripts/myNCL/NCLscripts/functions.ncl"

; Based on http://www.atmos.rcast.u-tokyo.ac.jp/nishii/programs/Plumb1985-grads.ncep.nov2014.txt
; Follows derivation of https://paperpile.com/view/3fd7d9e7-351f-08da-ba1d-af2f7191ed86 
; Calculated fluxes from climatological seasonal means
; with perturbations as differences from zonal mean
; Quasi-geostrophic, in spherical coordinates

begin

; Get experiment names

print("Calc Plumb fluxes")

lev_p = 250

a = 6.371e06	; radius of Earth
pi = 3.14159265358979

cp = 1.00464e3
Rd = 287.0
kappa = (Rd/cp)
omega =  7.2921e-5
g = 9.80665
P0 = 1000.0
dirin = "/esarchive/scratch/rwhite/data/Projects/seaice/ECEarth/"

dirout = dirin

filename = "SF250hPa_control_djf_ensmean.nc" 

file_in = addfile(dirin + filename,"r")

; check that latitudes are increasing
lats = file_in->lat
if (lats(0) .gt. lats(1) .or. lats(0) .gt. 0) then
    print(" latitude array is backwards from expected! Need to fix this")
    ;SF1 = file_in->SF(0,::-1,:)
    ;print(dimsizes(SF1))
    ;lats = file_in->lat(::-1)
    ; need to make sure it prints out correctly if we flip the inputs
else
    SF1 = file_in->SF(0,:,:)

end if

filo = "Plumb_" + filename
system("/bin/rm -f " + dirout + filo)
ncF = addfile(dirout + filo,"c")
fAtt = True            ; assign file attributes
fAtt@creation_date = systemfunc ("date")
fAtt@title = " Plumb fluxes calculated in NCL using printPlumb2D from script Calc_Plumb_flux_ECEarth.ncl"

returned = printPlumb2D(file_in,SF1,lev_p,ncF)
delete([/SF1,lats/])

end
